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1.  Introduction 


Direct- fire  projectiles  are  fired  by  line-of-sight  aiming  and  are  fired  from  ground-based 
platforms,  helicopters,  and  fixed  wing  aircraft.  A  number  of  conditions  can  cause  rounds  to  miss 
an  intended  target.  These  conditions  include  manufacturing  inaccuracies  of  the  gun  tube, 
propellant,  or  projectile,  along  with  variable  atmospheric  conditions,  firing  platfonn  motion,  and 
aiming  errors.  With  the  advent  of  low-cost,  small,  rugged,  microelectromechanical  systems, 
dramatic  reduction  of  dispersion  for  direct-fire  projectiles  equipped  with  a  relatively  inexpensive 
flight  control  system  is  possible.  One  design  concept  consists  of  a  set  of  controllable  canards 
located  near  the  nose  of  the  projectile.  This  report  develops  a  unique  flight  control  law  tailored 
to  the  control  of  smart  projectiles  through  the  application  of  model  predictive  control  and 
projectile  linear  theory. 

In  model  predictive  control,  a  dynamic  model  of  the  plant  is  used  to  project  the  state  into  the 
future  and  subsequently  use  the  estimated  future  state  to  determine  control  action.  It  has  been 
found  to  be  a  practical  and  increasingly  employed  control  technique  (1 ).  Currently,  model 
predictive  control  is  being  applied  to  a  wide  variety  of  problems,  spanning  many  different 
industries.  Mei  et  al.  (2)  studied  vibration  reduction  of  a  tall  building  experiencing  wind 
excitation  using  model  predictive  control  and  linear  quadratic  guassian  control  strategies.  They 
found  that  the  model  predictive  control  scheme  perfonned  well  and  was  robust  to  uncertainty  in 
building  stiffness.  Tsai  and  Huang  (3)  used  a  model  reference  adaptive  predictive  controller  for  a 
variable-frequency  oil-cooling  machine  used  in  concert  with  dynamically  complex  machine 
tools.  Kvatemik  et  al.  ( 4 )  developed  a  generalized  predictive  controller  for  tilt  rotor  aeroelastic 
stability  augmentation  in  airplane  mode  of  flight.  Using  the  model  predictive  control  strategy, 
significant  increases  in  damping  of  aircraft  body  vibration  modes  were  achieved  in  a  wind  tunnel 
test.  Siegers  and  Costello  (5)  applied  model  predictive  control  to  a  parafoil  for  autonomous 
delivery  of  a  payload  in  battlefield  conditions.  Burchett  and  Costello  ( 6)  used  a  simplified  form 
of  model  predictive  control  applied  to  a  projectile  with  lateral  pulsejets.  Their  strategy  was  to 
use  projectile  linear  theory  to  map  the  projected  impact  point  in  the  vertical  target  plane  and  base 
control  action  on  projected  miss  distance  and  direction.  The  key  difference  between  their 
strategy  and  that  detailed  here  is  that  the  control  strategy  employed  by  Burchett  and  Costello 
calculated  errors  only  in  the  target  plane,  while  the  control  strategy  used  here  considers  error 
along  the  length  of  the  trajectory.  In  addition,  the  pulsejet  control  scheme  used  by  Burchett  and 
Costello  is  inherently  discontinuous  and  can  only  be  applied  at  a  discrete  number  of  points. 
Canard  control,  as  applied  in  this  report,  is  continuous  and  is  applied  for  the  full  flight  duration. 

Any  model  predictive  control  scheme  is  dependent  upon  the  accuracy  of  the  underlying  dynamic 
model  representing  the  plant.  Under  most  flight  conditions,  the  equations  of  motion  for  a 
projectile  in  atmospheric  flight  can  be  adequately  represented  by  a  six-degree-of-freedom  rigid 
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body  model  with  externally-applied  aerodynamic  forces  and  moments.  The  resulting  differential 
equations  have  been  shown  to  provide  an  accurate  representation  of  projectile  flight 
characteristics  (7),  though  their  inherent  nonlinearity  prevents  direct  use  in  model-predictive 
control  applications.  However,  a  series  of  manipulations  and  simplifications  of  the  equations  of 
motion  allow  closed  form  solution  of  the  projectile  trajectory  under  restricted  flight  conditions. 
The  simplified  dynamic  equations  and  their  resulting  solutions  have  become  known  as  projectile 
linear  theory.  Projectile  linear  theory  has  been  extended  by  various  authors  to  handle  more 
sophisticated  aerodynamic  models  (5),  asymmetric  mass  properties  (9),  fluid  payloads  (10,  11), 
moving  internal  parts  (12, 13),  dual  spin  projectiles  (14, 15),  ascending  flight  (16),  and  lateral 
force  impulses  (1 7-20).  Aerodynamic  range  reduction  software  used  in  spark  range  facilities 
utilizes  projectile  linear  theory  in  estimation  of  aerodynamic  coefficients. 

The  reported  work  employs  model  predictive  control  and  projectile  linear  theory  for  control  of  a 
direct- fire  projectile.  The  basic  projectile  configuration  under  consideration  is  fin  stabilized,  and 
the  fins  are  slightly  canted  to  provide  moderate  roll  rates  during  flight.  A  set  of  controllable 
canards  located  near  the  nose  of  the  projectile  are  used  as  the  control  mechanism.  The  canards 
can  be  directed  to  provide  swerve  forces  and  pitch  and  yaw  moments  to  the  projectile.  The 
control  law  uses  an  approximate  closed  form  solution  of  projectile  motion  to  predict  the  states  of 
the  projectile  over  a  set  distance  known  as  the  prediction  horizon.  Current  and  future  control 
actions  are  determined  based  on  minimizing  the  estimated  error  of  future  states.  It  is  assumed 
that  sensor  feedback  is  provided  by  an  onboard  inertial  measurement  unit  (IMU).  Simulation 
results  to  establish  the  utility  of  the  new  model  predictive  flight  control  system  design 
methodology  are  generated  for  an  exemplar  projectile.  Parametric  trade  studies  are  conducted 
that  consider  the  effect  of  the  cost  function  weighting  matrices,  prediction  horizon  length,  state 
estimation  step  size,  and  the  model  update  interval  on  impact  point  dispersion. 


2.  Projectile  Dynamic  Model 


The  nonlinear  trajectory  simulation  used  in  this  study  is  a  standard  six-degree-of-freedom  model 
typically  used  in  flight  dynamic  modeling  of  projectiles.  A  schematic  of  the  projectile 
configuration  is  shown  in  figures  1  and  2.  The  six  degrees  of  freedom  are  the  three  inertial 
components  of  the  position  vector  from  an  inertial  frame  to  the  projectile  mass  center  and  the 
three  standard  Euler  orientation  angles. 

The  equations  of  motion  are  provided  in  equations  1-4  (21-23). 
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In  equations  1  and  2,  the  standard  shorthand  notation  for  trigonometric  functions  is  used: 
sin(«)  =  sa ,  cos(a)  =  ca  ,  andtan(a)  =  ta  .  The  forces  in  equation  3  contain  contributions  from 
weight  (W),  body  aerodynamics  (A),  and  the  control  canards  (C). 


Control  Canards 


Figure  1.  Schematic  of  the  position  coordinates  of  a  direct-fire  projectile. 
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The  dynamic  equations  are  expressed  in  a  body-fixed  reference  frame,  thus  all  forces  acting  on 
the  body  are  expressed  in  the  projectile  reference  frame.  The  weight  force  is  shown  inequation  6: 
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while  the  aerodynamic  force  acting  at  the  center  of  pressure  of  the  projectile  is  given  by  equation  7: 
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The  control  forces  are  the  aerodynamic  drag  forces  created  by  the  control  canards  in  the 
directions  perpendicular  to  the  axis  of  symmetry  of  the  projectile: 
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The  applied  moments  about  the  projectile  mass  center  contain  contributions  from  steady 
aerodynamics  (SA),  unsteady  aerodynamics  (UA),  and  the  control  canards  (C): 
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The  moment  components  due  to  steady  aerodynamic  forces  and  control  canard  forces  are 
computed  with  a  cross  product  between  the  distance  vector  from  the  mass  center  to  the  location 
of  the  specific  force  and  the  force  itself.  The  unsteady  body  aerodynamic  moment  provides  a 
damping  source  for  projectile  angular  motion  and  is  given  by  equation  10: 
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The  mass,  mass  center  location,  and  inertial  properties  of  the  projectile  are  all  assumed  to  be 
constant  throughout  the  duration  of  the  flight.  The  center  of  pressure  location  and  all 
aerodynamic  coefficients  (Cx0,  CYPA ,  CNA ,  CDD ,  CLP ,  CMQ  )  depend  on  local  Mach  number  and  are 

computed  during  simulation  using  linear  interpolation. 

The  dynamic  equations  given  by  equations  1—4  are  numerically  integrated  forward  in  time  using 
a  4th  order,  fixed  step  Runge-Kutta  algorithm.  Costello  and  Anderson  (7)  present  correlation  of 
this  dynamic  model  against  range  data  for  a  fin-stabilized  projectile. 
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3.  Projectile  Linear  Theory  Trajectory  Solution 


The  six-degree-of- freedom  rigid  body  projectile  model  shown  earlier  consists  of  12  highly 
nonlinear  differential  equations  for  which  a  closed  form  solution  has  not  been  directly  found. 
Significant  work  has  been  perfonned  to  simplify  the  equations  of  motion  such  that  an  accurate 
analytical  solution  can  be  determined.  In  order  to  arrive  at  a  set  of  analytically  solvable  ordinary 
linear  differential  equations,  the  following  assumptions  and  simplifications  are  made: 

1 .  Rather  than  employing  a  reference  frame  fixed  to  the  projectile  body,  projectile  linear 
theory  uses  an  intermediate  reference  frame  which  is  aligned  with  the  projectile  axis  of 
symmetry  but  does  not  roll.  Lateral  translational  and  rotational  velocity  components 
described  in  this  frame,  known  as  the  no-roll  frame  or  the  fixed  plane  frame,  are  denoted 
with  a  ~  superscript.  Components  of  the  linear  and  angular  body  velocities  in  the  fixed 
plane  frame  are  computed  from  body  frame  components  of  the  same  vector  through  a 
single  axis  rotation  transformation.  For  example,  the  body  frame  components  of  the 
projectile  mass  center  velocity  are  transformed  to  the  fixed  plane  by 

ii  10  0  u 

<  v  >=  0  c,  -s,  <  v  [.  (11) 


2.  A  change  of  variables  is  made  from  the  velocity  along  the  projectile  axis  of  symmetry,  u,  to 
the  total  velocity,  V.  Equations  12  and  13  relate  V  and  u  and  their  derivatives: 


V  =  \Ju^  +  \  +  w2  =\ju2  +  v2  +w2  ; 


uu  +  w  +  ww  uu  +  vv  +  ww 


3.  Dimensionless  arc  length,  s,  is  used  as  the  independent  variable  instead  of  time,  t. 
Equation  14  defines  dimensionless  arc  length: 

1  1 
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n  J 


Equations  15  and  16  relate  time  and  arc  length  derivatives  of  a  dummy  variable  C,.  Dotted  tenns 
refer  to  time  derivatives  and  primed  tenns  denote  arc  length  derivatives: 
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4.  Euler  pitch  and  yaw  angles  are  assumed  to  be  small  so  that 

sin(6*) «  6  cos(6*) « 1 

sin(^)  «  \ff  cos (y/)  ~  1 . 

5.  Aerodynamic  angles  of  attack  are  small  so  that 
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V  V 

6.  The  projectile  is  mass-balanced  such  that  the  center  of  gravity  lies  in  the  rotational  axis  of 
symmetry: 

4=4=4=° 

4=4  09) 

4=4=4  . 

7.  Quantities  l7  and  ^  are  large  compared  to  v’ vv’  and  r  such  that  products  of  small 
quantities  and  their  derivatives  are  negligible. 

A  more  detailed  discussion  of  the  development  of  projectile  linear  theory  is  provided  by  McCoy 
(21).  Application  of  the  above  stated  assumptions  leads  to  a  set  of  coupled  linear  differential 
equations,  with  the  exception  that  the  total  velocity,  V  ,  the  roll  rate,  p  ,  and  the  pitch  angle,  6 , 

appear  in  nonlinear  fashion  in  many  of  the  equations.  To  remedy  this,  the  assumption  is  made 
that  V  changes  slowly  with  respect  to  the  other  variables  and  is  thus  considered  to  be  constant, 

V  =  V0 ,  when  it  appears  as  a  coefficient  in  all  dynamic  equations  except  its  own.  In  addition,  the 
roll  rate  and  pitch  angle  are  held  constant,  p  =  p0  and  6  =  0O,  only  when  they  appear  in 
nonlinear  fashion.  The  equation  for  the  total  velocity  is  shown  as  equation  20: 
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The  remaining  1 1  equations  can  be  written  as 
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J 
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(31) 

(32) 


Closed  fonn  solutions  to  the  above  equations  can  be  found,  though  the  results  are  omitted  here. 

A  more  detailed  treatment  of  the  solutions  to  the  projectile  linear  theory  equations  can  be  found 
in  reference  (23). 

The  variables  CY0  and  Cz0  are  aerodynamic  trim  forces  perpendicular  to  the  projectile  axis  of 

symmetry,  which  are  created  by  movement  of  the  control  canards  and  are  treated  directly  as 
control  inputs. 

In  reality,  the  total  velocity,  V  ,  does  not  remain  constant  for  the  duration  of  the  flight;  therefore, 
the  total  velocity  must  be  periodically  measured  throughout  the  trajectory  and  updated  in  the 
remaining  equations.  The  center  of  pressure  location  and  the  aerodynamic  coefficients 
( Cx0 ,  CYPA ,  CNA ,  CDD ,  CLP ,  CM0  ),  which  all  depend  on  local  Mach  number,  must  also  be 

recomputed  each  time  V  is  updated.  The  effect  of  the  length  of  the  update  interval  on  the 
accuracy  of  the  model  was  studied  by  Burchett  et  al.  (19). 


4.  Model-Predictive  Flight  Control  System 


The  model-predictive  controller  uses  the  linearized  model  of  the  system  to  propagate  the  states 
forward  in  time  over  an  interval  known  as  the  prediction  horizon ^Hp  j  (24).  Control  action  is 

based  on  comparison  of  the  predicted  states  and  a  predetennined  desired  trajectory  over  the 
prediction  horizon.  As  the  prediction  step  is  marched  forward,  so  too  is  the  prediction  horizon — 
a  process  referred  to  as  the  “receding  horizon  principle.”  The  control  action  at  each  step  is 
determined  by  minimizing  a  quadratic  cost  function,  defined  as 

J  =  (W-f)T  Q(W-Y)  +  UTRU  .  (33) 

The  matrix  W  contains  the  desired  system  outputs,  w ,  over  the  length  of  the  prediction  horizon. 
The  desired  system  outputs,  w ,  at  each  prediction  step  consists  of  the  desired  x ,  y  ,  and  z 
coordinates  at  that  time  instant.  These  values  need  to  be  loaded  into  the  onboard  computer  prior 
to  projectile  launch. 
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In  order  to  develop  an  expression  for  the  predicted  system  outputs  over  the  prediction  horizon, 
the  system  is  first  cast  in  standard  discrete  state-space  fonn: 

xk+l  =  A(As)xk  +  B(A5,)w/t  +  F(A5)  (36) 

yk=Cxk, 

where  the  values  within  the  matrices  A,  B  ,  and  F  depend  on  the  arc  length  step  size(  A  s).  The 
projectile  linear  theory  expressions  shown  in  the  previous  section  are  used  to  form  the  state  space 
matrices  through  a  14-step  loop  in  the  control  algorithm.  In  the  first  step,  all  the  states  and 
controls  are  set  to  zero  and  the  solutions  are  evaluated  over  one  arc  length  step  to  detennine  the 
values  within  the  constant  vector  F  .  In  the  next  step,  the  first  state,  u  ,  is  set  equal  to  one,  with 
the  remaining  states  and  controls  still  equal  to  zero,  and  the  expressions  are  reevaluated.  By 
subtracting  the  values  of  the  constants,  F  ,  the  coefficients  making  up  the  first  column  of  A  can 
be  found.  This  process,  consisting  of  setting  a  state  variable  equal  to  one,  evaluating  the  linear 
theory  solutions,  then  subtracting  the  constant  values,  is  repeated  for  each  of  the  remaining  10 
states  to  fully  populate  the  state  matrix  A  one  column  at  a  time.  The  control  matrix,  B  ,  is 
formed  in  exactly  the  same  manner  with  all  1 1  states  equal  to  zero  and  the  controls,  CY0  and  Cz0 
alternately  set  equal  to  one. 


The  desired  outputs  of  the  system  are  its  center  of  mass  position  states  (x,  y,  and  z).  The  matrix 
C  is  then  simply 

"0  000100000  0~ 

C=00000  1  00000.  (37) 

00000010000 
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A  recursive  formula  can  be  found  for  yk+j,  1  <  j  <  H p  by  substituting  the  expression  for  xk+j  into 
the  expression  for  yk+j .  The  result  is 

=CA+Se^'%.,, +'Zca,-‘f,  (38) 

i= 1  i= 1 


or,  in  matrix  form. 


Y  —  KCAxk  +  KcabU  +  Kcaf  , 


(39) 


where 
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0  ...  0 

o  •••  ; 

CB  ■■■  0 

;  ■•.  o 
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,  and 
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CAP 


CF 

CAF  +  CF 
CA2F  +  CAF  +  CF 


CA  HP~'F  +  ---  +  CAF  +  CF 


(40) 


(41) 


(42) 


The  cost  function,  J ,  then  becomes 

J  =  (W- KCAxk  - KcabU  - Kcaf)t Q(W  - KCAxk  - KcabU -Kcaf)  +  UtRU.  (43) 

The  minimum  of  the  cost  function  is  determined  by  selecting  the  control  input  vector  that  forces 
the  gradient  of  the  cost  function  to  zero: 

U  =  K(W-Kaxt-KCJP),  (44) 

where 

K=(Kca/QKca,+R)'KcCQ-  <45> 
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It  should  be  noted  that  U  contains  the  optimal  control  inputs  over  the  entire  prediction  horizon. 
At  each  arc  length  step,  k ,  only  uk  is  used,  which  is  the  first  element  of  U  .  The  first  element  of 

U  is 


u,=K,(W-Kax,-Kar),  (46) 

where  Kl  consists  only  of  the  first  M  rows  of  K.  Note  that  M  is  defined  as  the  number  of 
control  inputs  which,  in  this  application,  is  two  ( CY0  and  Czo ). 

It  is  assumed  that  full-state  feedback  is  available  for  use  in  the  control  law;  that  is, 
x  ,y,z,i//,0,c/>,u,v,w,p,q,  and  r  are  sensed  or  estimated  by  the  IMU.  Furthermore,  the  weapon 
that  fires  the  projectile  provides  a  desired  trajectory  leading  to  the  target.  At  time  =  0,  the 
controller  is  provided  with  the  full  state  of  the  projectile.  The  total  velocity,  V  ,  is  calculated 
from  the  projectile  mass  center  velocity  states  and  set  to  V0  in  the  linear  model.  The  linear 
model  is  then  used  to  propagate  the  remaining  1 1  states  forward  by  A  s.  These  values  are  used  to 
populate  the  A,  B  ,  and  A  matrices,  which  are  sent  to  the  MPC  routine.  The  MPC  routine 
calculates  the  optimal  control  sequence  over  the  length  of  the  update  interval.  When  the 
projectile  has  covered  the  length  of  the  first  update  interval,  as  well  as  every  subsequent  update 
interval,  the  controller  is  provided  with  full  state  feedback  and  the  process  is  repeated.  The 
control  sequence  calculated  by  the  model  predictive  controller  contains  control  inputs  at 
increments  of  as.  Linear  interpolation  is  applied  to  determine  control  inputs  between 
increments  of  A  s. 

It  is  important  to  note  that  the  controls  resulting  from  the  earlier  calculations  are  expressed  in  the 
fixed  plane  frame,  as  per  assumption  number  1  in  the  linear  theory  section  of  this  report.  To  be 
applied  to  the  canards,  the  control  inputs  must  be  converted  to  the  conventional  body  fixed 
reference  frame. 

Each  time  the  Mach  number  is  updated  in  the  linear  model,  the  matrices  A,  B,  and  F  are  updated 
as  well.  This  in  turn  requires  updating  of  the  gain  matrices.  The  size  of  each  of  these  matrices, 
and  hence  the  computational  time  required  to  calculate  them,  is  governed  by  the  length  of  the 
prediction  horizon.  Obviously,  frequent  updates  to  the  linear  model  and  a  long  prediction 
horizon  provide  greater  accuracy  in  the  predictor  and  more  efficient  control.  These  observations 
are  tempered  with  the  need  to  limit  the  computational  demand  placed  on  the  onboard  processor. 


5.  Results 


To  establish  the  utility  of  the  model-predictive  controller  in  a  projectile  application,  a  4.5-ft  long 
fin-stabilized  projectile  is  considered.  The  projectile  has  a  total  weight  of  22.9748  lb,  a  center  of 
gravity  location  of  2.5  ft  from  the  base,  and  four  tail-mounted  stabilization  fins.  The  roll  and 
pitch  inertia  of  the  body  are  0.0057  and  1.35  slug-ft2,  respectively.  A  set  of  controllable  canards, 
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which  alter  the  aerodynamic  forces  and  moments,  are  located  4.25  ft  from  the  base  of  the 
projectile. 

To  model  uncertainty  in  launch  conditions,  which  is  a  primary  cause  of  dispersion,  the  initial 
pitch  and  yaw  rates,  pitch  and  yaw  angles,  and  body  velocities  are  all  considered  to  be  normally- 
distributed  random  numbers  with  means  and  standard  deviations  that  are  representative  of  actual 
launch  uncertainties.  The  values  chosen  are  shown  in  table  1 . 


Table  1.  Initial  condition  uncertainty  parameters  for  dispersion  analysis. 


Initial  Condition 

Mean 

Standard  Deviation 

Pitch  angle  (9) 

0.1  rad 

0.01  rad 

Yaw  angle  (v|/) 

0.0  rad 

0.01  rad 

Pitch  rate  (q) 

-0.18  rad/s 

2.0  rad/s 

Yaw  rate  (r) 

0.0  rad/s 

2.0  rad/s 

x  body  velocity  (u) 

1143.3797  ft/s 

15  ft/s 

y  body  velocity  (v) 

-0.00002502  ft/s 

3  ft/s 

z  body  velocity  (w) 

0.375346  ft/s 

3  ft/s 

The  desired  trajectory  is  chosen  as  that  which  the  projectile  would  follow  in  the  absence  of 
uncertainty  with  initial  conditions  of  x0,  y0,  z0,  y/0,  <f>0,  and  r0  =0,  00=  0.1  rad,  u0  =  1 143.38 
ft/s,  v0  =  -2.502x1  O'5  ft/s,  w0=  0.375  ft/s,  p0  =  51.5  rad/s,  and  q0  =  -0.18  rad/s.  The  target 
location  is  chosen  as  x=  6216.613  ft,  y  =  0.261  ft,  and  z  =  0.0  ft.  Figure  3  shows  typical 
dispersion  results  for  50  sample  trajectories  with  no  control  applied  and  initial  condition 
perturbations  as  described  earlier.  The  circular  error  probable  (CEP)  shown  in  the  figure  is  based 
on  a  50%  hit  criterion;  that  is,  the  CEP  is  defined  as  the  minimum  radius  of  a  circle  centered  at 
the  mean  impact  point  and  containing  at  least  50%  of  the  shot  impact  points.  With  no  control 
applied,  the  CEP  is  106  ft.  For  reference,  a  second  CEP  is  shown,  which,  instead  of  being 
centered  at  the  mean  impact  point,  is  instead  centered  at  the  target  location.  A  50%  hit  criterion 
is  still  used.  The  second  CEP  has  a  radius  of  1 13.5  ft. 


Figure  4  shows  the  dispersion  results  with  model-predictive  control  applied.  The  prediction 
horizon,  H  ,  is  chosen  as  50.  The  error  weighting  matrix,  Q ,  is  chosen  to  be  a  function  of  range 

as  follows: 


Uoooj 


Q  =  q 


0 

0 


0 


Uoooj 

0 


0  0 

0  0 


Jk+H„- 1 


1000 


(47) 
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Figure  3.  Uncontrolled  dispersion  (CEP  =  113.5  ft  centered  at  mean  impact  point). 

in  which  q  is  chosen,  in  this  case,  to  be  0.5.  By  defining  Q  in  this  manner,  the  error  weighting 
is  increased  quadratically  as  the  projectile  flies  downrange.  This  prevents  the  tendency  for  the 
controller  to  attempt  to  force  the  projectile  onto  the  desired  trajectory  immediately  after  launch, 
leading  to  a  large  initial  control  input  followed  by  subsequent  control  inputs  of  nearly  zero.  In 


addition,  the  control  weighting  matrix, 

R, 

is  defined 

to  be 

~hp 

0 

0 

0 

0" 

0 

Hp~1 

0 

0 

0 

R  =  r 

0 

0 

0 

0 

2 

0 

0 

0 

0 

1 

where  r  is  chosen  as  2.0  in  this  case.  By  defining  R  in  this  manner,  the  current  control  value  is 
weighted  H  times  heavier  than  the  control  value  at  the  end  of  the  prediction  horizon.  This 
prevents  large  controls  from  being  chosen  at  the  beginning  of  an  update  interval,  even  if 
significant  error  is  present.  The  model  update  interval  for  the  case  shown  in  figure  2  is  1000  arc 
lengths,  and  the  arc  length  step  size,  a  s,  is  20.  These  parameters  provide  a  good  baseline  from 
which  to  begin  examining  the  performance  of  the  model-predictive  controller.  Figures  5  and  6 
show  a  typical  controlled  and  uncontrolled  trajectory  with  the  model  predictive  control 
parameters  set  listed  earlier. 
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Figure  4.  Controlled  dispersion  (CEP  =  0.02  ft  centered  at  mean  impact  point). 


Figure  5.  Typical  altitude  vs.  range. 
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Figure  6.  Typical  cross  range  vs.  range. 

The  model-predictive  controller  provides  a  very  significant  reduction  in  the  CEP — from  106  ft, 
in  the  uncontrolled  case,  to  0.02  ft  (or  less  than  1/4  in)  with  control  applied.  It  should  also  be 
noted  that  the  mean  impact  point  is  almost  0.06  ft,  or  -0.75  in,  above  the  target  location.  In  the 
cross-range  direction,  however,  the  mean  impact  point  is  only  0.0085  ft  away  from  the  target. 
This  bias  error  in  the  z-direction  can  be  attributed  directly  to  errors  in  the  linear  model  used  in 
the  predictor.  One  of  the  primary  assumptions  upon  which  projectile  linear  theory  is  based  is 
that  the  projectile  maintains  a  small  angle  of  attack.  As  the  target  is  approached,  the  angle  of 
attack  of  the  projectile  is  forced  to  a  small,  nonzero  number.  Though  it  isn’t  necessarily  in 
violation  of  the  small  angle  of  attack  assumption,  it  is  enough  to  cause  a  small  deviation  between 
the  trajectory  predicted  with  the  linear  model  and  that  which  is  arrived  at  by  integrating  the  full 
six-degree-of-freedom,  nonlinear  equations.  This  error  is  demonstrated  by  plotting  the  error 
between  the  validated,  full  six-degree-of-freedom,  nonlinear  trajectory,  which  is  solved  using  a 
fixed  step,  4th  order  Runge-Kutta  method,  and  the  linear  theory  trajectory  solution.  The  linear 
solution  is  corrected  to  match  the  nonlinear  solution  every  1000  arc  lengths  to  mimic  flight 
control  system  feedback.  The  control  input  is  set  equal  to  zero  in  both  cases  and  the  initial 
conditions  are  set  to  match  those  used  in  creating  the  desired  trajectory.  Figure  7  shows  the 
linear  theory  error  as  a  function  of  arc  length  in  the  x,  y,  and  z  directions. 

Note  that  the  error  is  of  the  same  order  of  magnitude  near  the  end  of  the  trajectory  as  the 
variation  of  the  mean  impact  point  in  the  CEP  plot.  It  should  also  be  noted  that  the  error  is 
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Figure  7.  Error  between  linear  and  nonlinear  trajectory  solutions. 

greatest  at  the  beginning  and  end  of  the  flight,  where  the  trajectory  is  furthest  from  horizontal. 

At  the  midpoint  of  the  trajectory,  where  the  path  of  flight  is  nearly  flat  and  the  projectile  angle  of 
attack  is  nearly  zero,  the  error  in  all  three  spatial  directions  also  becomes  very  close  to  zero. 

Figure  8  shows  the  required  control  inputs  for  the  trajectory  shown  in  figures  3  and  4.  The 
magnitudes  of  the  control  inputs  required  to  achieve  the  shown  degree  of  tracking  are  attainable 
for  a  set  of  nose-mounted  canards. 

Great  care  must  be  taken  when  choosing  the  gain  values,  q  and  r ,  such  that  the  control 
values,  CY0  and  Czo ,  never  exceed  approximately  one.  Such  large  control  inputs  violate  the 
small  angle  of  attack  assumption  upon  which  the  linearized  model  is  based.  As  a  result,  the 
linear  model  no  longer  accurately  approximates  the  true,  nonlinear  system,  and  the  controller 
loses  its  ability  to  accurately  predict  future  states.  When  the  system  is  provided  with  state 
feedback  from  the  IMU  under  these  circumstances,  there  are  very  large  errors  and  the  controller 
subsequently  attempts  to  choose  a  large  control  value  to  compensate.  Within  one  to  two  update 
intervals,  the  error  becomes  large  enough  that  control  saturates. 

Other  applications  of  model  predictive  control,  such  as  that  discussed  by  Mei  et  al.  (2),  use  an 
iterative  scheme  to  set  a  maximum  control  input  value.  In  the  application  being  discussed,  where 
speed  of  control  computation  is  extremely  critical,  iteration  is  not  practical.  If  the  processor  is 
occupied  by  an  iterative  routine  while  the  projectile  continues  to  fly  downrange,  control  is  lost. 
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Figure  8.  Required  control  inputs  for  a  typical  trajectory. 

A  second  option  is  to  simply  clip  the  control  values  at  the  maximum  allowable  value.  This,  too, 
presents  problems  as  future  controls  are  calculated  under  the  assumption  that  all  previous 
controls  were  applied  exactly  as  calculated.  When  this  clipping  scheme  is  attempted,  the  control 
begins  oscillating  rapidly  between  both  allowable  extremes,  quickly  resulting  in  instability. 

The  results  discussed  earlier  assume  perfect  sensor  feedback,  which  in  reality  can  never  be 
achieved.  In  practice,  sensors  possess  error  created  by  both  bias  and  noise.  These  are  modeled 
in  the  simulation  by  choosing  normally-distributed  random  numbers  with  means  of  zero  and 
standard  deviations  that  are  representative  of  commonly  used  IMU  sensors.  A  bias  value  is 
randomly  chosen  for  each  sensor  at  the  start  of  every  flight  simulation  and  retained  throughout 
that  particular  simulation.  In  addition,  a  noise  value  is  randomly  chosen  for  each  sensor  every 
time  feedback  is  implemented.  Both  the  bias  and  noise  value  are  added  to  the  sensor  readings  at 
each  update  interval.  The  sensor  bias  and  noise  standard  deviations  used  are  summarized  in 
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table  2.  All  subsequent  results  displayed  in  this  report  employ  sensor  bias  and  noise  applied  in 
this  way. 

Table  2.  Sensor  noise  and  bias  values. 


Sensor  Function 

Bias  Standard  Deviation 

Noise  Standard  Deviation 

x  and  y  position  (ft) 

0.52 

0.52 

z  position  (ft) 

1.18 

1.18 

x  and  y  velocity  (ft/s) 

0.10 

0.08 

z  velocity  (ft/s) 

0.16 

0.13 

roll,  pitch,  and  yaw  angles  (deg) 

0.3 

0.3 

roll,  pitch,  and  yaw  rates  (deg/s) 

0.05 

0.01 

Sensor  noise  and  bias  become  the  dominant  sources  of  error  when  they  are  applied  in  this 
application.  As  the  standard  deviation  of  the  sensor  noise  remains  constant  throughout  the 
projectile  flight,  it  no  longer  makes  sense  to  define  the  error  weighting  matrix  as  a  function  of 
projectile  range.  Doing  so,  while  keeping  the  value  of  q  low  enough  to  avoid  violating  the  small 
angle  of  attack  assumption  near  the  end  of  the  trajectory,  unnecessarily  limits  the  control  action 
near  the  beginning  of  the  trajectory.  The  error  weighting  matrix,  Q ,  is  instead  defined  simply  as 
the  identity  matrix  multiplied  by  the  constant  gain  value,  q  . 

Figure  9  shows  dispersion  results  with  sensor  noise  and  bias  applied.  The  prediction  horizon,  Hp, 
is  again  50.  The  error  gain,  q,  is  set  to  1.  The  control  weighting  matrix,  R ,  is  defined  as  shown 
in  equation  48  with  the  gain  value,  r ,  equal  to  0.2.  The  update  interval  is  1000  arc  lengths  and 
the  arc  length  step  size,  a  s,  is  20.  As  a  direct  result  of  the  sensor  uncertainty,  the  CEP  radius  is 
increased  to  2.1  ft. 

In  figure  10,  the  effects  of  changing  both  the  control  gain,  r ,  and  the  prediction  horizon,  Hp,  are 
shown.  The  model  update  interval  is  held  constant  at  1000  arc  lengths  and  the  arc  length  step 
size  is  held  constant  at  20  arc  lengths  throughout  all  of  the  simulations  shown  in  figure  10.  In 
addition,  the  error  weighting  matrix,  Q,  is  defined  as  the  identity  matrix  with  the  error  gain,  q, 
equal  to  1.  The  control  weighting  matrix,  R ,  is  again  defined  as  shown  in  equation  48.  The 
control  gain,  r,  is  varied  from  0.025  to  10  with  the  prediction  horizon  held  constant.  The  process 
is  repeated  four  times  with  values  of  H p  =  25 ,  H p  =  50  ,  H p  =  75  ,  and  H p  =  100 . 

As  the  value  of  r  is  increased,  additional  weight  is  given  to  the  value  of  the  control  in  the  cost 
function  (equation  43).  This  in  turn  forces  the  magnitude  of  the  chosen  control  values  to  be 
smaller,  which  provides  less  control  authority.  As  would  be  expected,  figure  10  shows  that 
larger  values  of  r  lead  to  increased  dispersion.  However,  there  is  a  value  of  r  below  which  the 
control  values  are  allowed  to  be  too  large,  leading  to  violation  of  the  small  angle  of  attack 
assumption  and  loss  of  control.  This  minimum  value  of  r  varies  depending  on  the  length  of  the 
prediction  horizon.  In  figure  10,  the  lowest  attempted  values  of  r  which  resulted  in  a  controllable 
trajectory  are  shown  as  the  first  data  point  for  each  series. 
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Figure  9.  Controlled  dispersion  with  sensor  bias  and  noise  applied  (CEP  =  2.1  ft 
centered  at  mean  impact  point). 


Figure  10.  Controlled  dispersion  results  as  the  control  gain,  r,  is  varied. 
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It  is  also  apparent  from  figure  10  that,  for  a  given  value  of  r,  there  is  a  direct  relationship 
between  the  length  of  the  prediction  horizon  and  the  amount  of  impact  point  dispersion. 

Allowing  the  controller  to  take  into  account  an  increased  number  of  the  predicted  states,  as  a 
longer  prediction  horizon  does,  leads  to  more  intelligent  control  choices.  It  also  significantly 
increases  the  amount  of  computation  required  at  each  update  interval,  necessitating  a  more 
expensive  onboard  processor. 

A  similar  study  was  performed  to  investigate  the  effect  of  the  length  of  the  linear  model  update 
interval  on  the  dispersion  radius.  The  error  weighting  matrix,  Q,  is  again  the  identity  matrix  with 
the  error  gain,  q,  set  to  one.  The  control  weighting  matrix,  R,  is  defined  as  shown  in  equation  48 
with  the  gain  value,  r,  equal  to  0.5.  The  arc  length  step  size  is  again  20  arc  lengths.  The  linear 
model  update  interval  is  varied  from  100  to  2000  while  holding  the  value  of  the  prediction 
horizon  constant.  As  before,  the  process  is  repeated  four  times  with  prediction  horizon  values  of 
H p  =  25 ,  H p  =  50  ,  Hp  =  75  ,  and  Hp  =  100  .  The  results  are  shown  in  figure  11. 


Figure  11.  Controlled  dispersion  results  as  the  linear  model  update 
interval  is  varied. 

Longer  linear  model  update  intervals  lead  to  an  increase  in  dispersion.  This  increase  becomes 
more  apparent  for  update  intervals  greater  than  1000  arc  lengths.  For  prediction  horizon  lengths 
of  25  and  50  arc  length  steps,  update  intervals  greater  than  1600  arc  lengths  led  to  a  loss  of 
control.  This  results  from  the  linearized  model  deviating  too  far  from  the  true,  nonlinear  system. 
Upon  update,  the  error  becomes  very  large  and  a  large  control  is  chosen  to  compensate,  which  in 
turn  violates  the  small  angle  of  attack  assumption  and  causes  further  error  in  the  linear  model. 
Reducing  the  error  gain,  q,  or  increasing  the  control  gain,  r,  would  prevent  this  scenario  from 
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occurring.  However,  the  tradeoff  would  be  a  reduction  in  control  authority  and  an  increase  in 
dispersion. 

The  final  study  investigates  the  effects  of  the  length  of  the  step  size,  a  s,  used  by  the  controller  to 
propagate  the  linear  model  forward.  The  prediction  horizon  is  held  constant  at  50  steps  and  the 
error  weighting  matrix  is  the  identity  matrix  with  the  gain,  q  ,  equal  to  one.  The  control 
weighting  matrix  is  defined  as  shown  in  equation  48.  Four  values  for  the  arc  length  step  size 
( A s  =  5,  A s  =  10,  A  s  =  20  ,  and  A  s  =  40)  are  used  while  the  control  gain,  r ,  is  varied  over  the 
range  which  provided  suitable  control  inputs  for  each  step  size.  The  results  can  be  seen  in 
figure  12. 


Figure  12.  Controlled  dispersion  results  as  the  arc  length  step  size  is 
varied. 

As  with  any  discrete,  linear  model,  the  length  of  the  step  size  has  no  effect  on  the  accuracy  of  the 
model  itself.  However,  as  evidenced  by  figure  12,  the  length  of  the  arc  length  step  size  does 
have  an  effect  on  the  overall  accuracy  of  the  controller.  This  results  from  an  interplay  between 
two  competing  effects.  The  prediction  horizon  length  is  measured  in  the  number  of  steps  into  the 
future  that  are  used  in  the  calculation  of  the  optimal  control.  Therefore,  for  a  given  prediction 
horizon  length,  increasing  the  arc  length  step  size  allows  the  predictor  to  take  into  account  state 
values  farther  into  the  future.  However,  control  values  are  only  calculated  at  each  step 
increment,  with  control  values  between  calculation  steps  derived  from  linear  interpolation.  A 
large  step  size  can  therefore  lead  to  a  decrease  in  resolution  of  the  controller.  For  arc  length  step 
sizes  of  5,  10,  and  20,  these  effects  do  little  more  than  change  the  acceptable  range  of  gains, 
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shifting  the  lines  to  the  right  on  figure  12  for  increasing  values  of  A  s.  However,  at  A  s  =  40,  the 
controller  becomes  unable  to  provide  the  necessary  amount  of  oscillation  and  dispersion  is 
increased  dramatically. 


6.  Conclusions 


This  report  develops  a  method  for  applying  model  predictive  control,  a  proven  and  effective 
control  technique,  to  a  smart  projectile  application.  The  control  law  is  shown  to  dramatically 
reduce  the  impact  point  dispersion  caused  by  launch  disturbances.  The  method  uses  full  state 
feedback  to  create  a  linearized  model  of  the  projectile  and  quickly  predict  the  future  states  of  the 
system.  These  calculations  can  be  performed  by  a  relatively  inexpensive  onboard  processor.  As 
the  predicted  states  depend  on  the  states  provided  by  the  feedback  loop,  sensor  accuracy  is  very 
important  to  the  performance  of  the  system  and  was  shown  here  to  be  the  limiting  factor  in 
dispersion  reduction. 

Considerable  opportunities  exist  for  the  control  system  designer  to  tune  the  model-predictive 
controller  based  on  the  desired  application.  It  was  shown  that  the  length  of  the  prediction 
horizon  has  a  considerable  effect  on  the  dispersion  radius,  with  a  longer  prediction  horizon 
leading  to  a  decrease  in  dispersion.  However,  a  longer  prediction  horizon  increases  the  size  of 
the  matrices  used  in  the  control  calculation,  which  subsequently  necessitates  an  increase  in  the 
processing  power  required  to  perfonn  control  calculations  in  a  sufficiently  short  period  of  time. 
Shorter  linear  model  update  intervals  lead  to  a  decrease  in  dispersion  as  well,  but  with  a  similar 
increase  in  the  amount  of  onboard  computation  required.  The  length  of  the  arc  length  step  size 
was  shown  to  have  little  effect  on  dispersion  as  long  as  it  remained  below  20  arc  lengths. 

Control  and  error  gains  should  be  adjusted  to  allow  sufficient  control  authority  without  violating 
any  of  the  assumptions  upon  which  linear  theory  is  based.  No  iterative  scheme  is  built  into  the 
controller  to  limit  the  size  of  the  control  inputs,  so  the  control  system  designer  should  run  a  series 
of  simulations  prior  to  launching  a  projectile  to  ensure  that  the  control  and  error  gains  are 
properly  adjusted. 
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List  of  Symbols,  Abbreviations,  and  Acronyms 


P: 


r  ■ 

^  NA  * 


C 


MQ  • 


c  ■ 


C  C  ■ 

WO’  '“'ZO  ' 


C  : 


C  : 

^YPA 


C  : 

LP 


D  ■ 


Ip’^R’ 


m  : 

p,  q,  r  : 

ASL  '■ 
ASL  : 

m 

A  SLC: 
u,  v,  w  : 


Air  density 

Normal  force  aerodynamic  coefficient 

Pitch  rate  damping  moment  aerodynamic  coefficient 

Aerodynamic  drag  coefficient  in  direction  parallel  to  projectile  motion 

Aerodynamic  trim  coefficients  perpendicular  to  projectile  axis  of  symmetry 

Roll  moment  from  fin  cant 

Magnus  force 

Roll  damping  moment 

Projectile  characteristic  diameter 

Projectile  precesional  and  rotational  inertia 

Projectile  mass 

Angular  velocity  vector  components  expressed  in  the  body  fixed  reference 
frame 

Euler  yaw,  pitch,  and  roll  angles 

Stationline  distance  from  the  projectile  center  of  pressure  location  to  the  CG 

Stationline  distance  from  the  projectile  Magnus  force  location  to  the  CG 

Stationline  distance  from  the  control  canard  location  to  the  CG 

Translation  velocity  components  of  the  projectile  center  of  mass  resolved  in  the 
body  fixed  reference  frame 
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X,  y,  z  : 

V  ■ 

L,M,N : 

X,Y,Z: 

Hp'- 
s : 


Position  vector  components  of  the  projectile  mass  center  expressed  in  the 
inertial  reference  frame 

Magnitude  of  the  mass  center  velocity 

Total  external  applied  moment  on  the  rocket  about  the  mass  center  expressed  in 
the  rocket  reference  frame 

Total  external  applied  force  on  the  rocket  expressed  in  the  rocket  reference 
frame 

Prediction  horizon  used  in  model  predictive  controller 
Dimensionless  arc  length 
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PANAMA  CITY  FL  32407-7001 


2  CDR  NSWC 
J  FRAYSEE 
D  HAGEN 

17320  DAHLGREN  RD 
DAHLGREN  VA  22448-5000 

5  CDR  NSWC 
INDIAN  HEAD  DIV 

D  GARVICK  CODE  40D 
L  FAN  CODE  41 10C 
V  CARLSON  CODE  4120 
HLAST  CODE  4140E 
T  GRIFFIN  CODE  45  0D 
101  STRAUSS  AVE 
INDIAN  HEAD  MD  20640-5000 

1  CDR  NSWC 
INDIAN  HEAD  DIV 
LIBRARY  CODE  8530 
BLDG  299 

101  STRAUSS  AVE 
INDIAN  HEAD  MD  20640 

2  US  MILITARY  ACADEMY 
MATH  SCI  CTR  OF  EXCELLENCE 
DEPT  OF  MATHEMATICAL  SCI 
MDN  A 

MAJ  D  ENGEN 
R  MARCHAND 
THAYER  HALL 
WEST  POINT  NY  10996-1786 

3  CDR  US  ARMY  YUMA  PG 
STEYP  MT  AT  A 

A  HOOPER 
STEYP  MT  EA 
YUMA  AZ  85365-9110 

6  CDR  NSWC 
INDIAN  HEAD  DIV 
CODE  570D  J  BOKSER 
CODE  5710  L  EAGLES 
J  FERSUSON 

CODE  57  C  PARIS 

CODE  5710G  SKIM 

CODE  5710E  S  JAGO 

101  STRAUSS  AVE  ELY  BLDG 

INDIAN  HEAD  MD  20640-5035 

1  MICHIGAN  STATE  UNIV 
B  KIM 

2120  ENGINEERING  BLDG 
EAST  LANSING  MI  48824-1226 
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2  INDUSTRIAL  OPERATION  CMD 
AMFIO  PM  RO 

W  MCKELVIN 

MAJ  BATEMAN 

ROCK  ISLAND  IL  61299-6000 

3  PROGRAM  EXECUTIVE  OFFICER 
TACTICAL  AIRCRAFT  PROGRAMS 
PMA  242  1 

MAJ  KIRBY  R242 
PMA  242  33 
R  KEISER  (2  CPS) 

1421  JEFFERSON  DAVIS  HWY 
ARLINGTON  VA  22243-1276 

3  ARROW  TECH  AS  SOCIATES  INC 

R  WHYTE 
A  HATHAWAY 
H  STEINHOFF 

1233  SHELBOURNE  RD  STE  D8 
SOUTH  BURLINGTON  VT  05403 

3  US  ARMY  AVIATION  CTR 
DIR  OF  COMBAT  DEV 
ATZQ  CDM  C 
B  NELSON 
ATZQ  CDC C 
T  HUNDLEY 
ATZQ  CD 
G  HARRISON 
FORT  RUCKER  AL  36362 

1  US  ARMY  CORPS  OF  ENGINEERS 
CEERD  ZA  T 
W  MORRISON 
7701  TELEGRAPH  RD 
KINGMAN  BLDG  RM  322 
ALEXANDRIA  VA  22315 

1  DIRCTRT  OF  TRAINING,  DOCTRINE, 
AND  COMBAT  DEV 
ATZK  TDD  ORSA 
A  POMEY 

RM  202  HARRIS  HALL 

BLDG  2368  OLD  IRONSIDES  AVE 

FORT  KNOX  KY  40121-5000 


ABERDEEN  PROVING  GROUND 

3  CDR 

USA  ARDEC 
AMSTA  AR  FSF  T 
R  LIESKE 
J  WHITESIDE 
J  MATTS 
BLDG  120 

1  CDR 

USA  TECOM 
AMSTE  CT 
T  J  SCHNELL 
RYAN  BLDG 

3  CDR 

USA  AMSAA 
AMXSY  EV 
G  CASTLEBURY 
R  MIRABELLE 
AMXSY  EF 
S  MCKEY 

43  DIR  USARL 

AMSRD  ARL  WM 
T  ROSENBERGER 
AMSRD  ARL  WM  BA 
W  HORST  JR 
W  CIEPELLA 
F  BRANDON 
T  BROWN  (5  CPS) 

L  BURKE 
J  CONDON 
B  DAVIS 

T  HARKINS  (5  CPS) 

D  HEPNER 
V  LEITZKE 
M  HOLLIS 
A  THOMPSON 
G  BROWN 

AMSRD  ARL  WM  BB 
B  HAUG 

AMSRD  ARL  WM  BC 
P  PLOSTINS  (4  CPS) 

G  COOPER 
B  GUIDOS 
JSAHU 
M BUNDY 
D  LYON 

AMSRD  ARL  WM  BC 
J  BENDER 
J  GARNER 
S  WILKERSON 
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W  DRYSDALE 
R  COATES 
A  MIKHAL 
J  WALL 

AMSRD  ARL  WM  BD 
B  FORCH 

AMSRD  ARL  WM  BE 
M  SCHMIDT 
AMSRD  ARL  WM  BF 
AMSRD  ARL  WM  BR 
C  SHOEMAKER 
J  BORNSTEIN 
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